Plotting function to reduce code length + simplify. Then, short function to relabel MLS weight smoothing kernel to make visually consistent, updates are:

  • 3.0 –> 3.6
  • 4.0 –> 4.8
  • 5.0 –> 6.0
  • 6.0 –> 7.2
  • 7.0 –> 8.4
create_plot <- function(dataset, y_est, type, thresh, y_min, y_max, y_lab, add_title = TRUE) {
  plot <- dataset %>%
    ggplot(aes(x = .data[[type]], y = .data[[y_est]], fill = .data[[type]], color = .data[[type]])) +
    geom_rain(alpha = .5, rain.side = 'l',
              boxplot.args = list(color = "black", outlier.shape = NA),
              boxplot.args.pos = list(
                position = ggpp::position_dodgenudge(x = .1), width = 0.1
              )) +
    facet_grid(~study) +
    ylab(y_lab) +
    ylim(y_min,y_max)+
    theme_classic() 
  
  if (add_title) {
    plot <- plot + ggtitle(paste("Distribution by", type, "category (", thresh, "mask)"))
  }
  
  plot <- plot +
    scale_fill_manual(values = pal) +
    scale_color_manual(values = pal) +
    guides(fill = 'none', color = 'none') +
    theme(text = element_text(family = "Times New Roman"),
          axis.text = element_text(size = 12, angle = 45, hjust = 1),
          axis.title = element_text(size = 12),
          legend.text = element_text(size = 12),
          legend.title = element_text(size = 12),
          plot.title = element_text(size = 16))
}

by_sample_conmod <- function(dataset, y_est, y_min, y_max, y_lab) {
  
  # Define a function to create individual plots
  create_plot <- function(data, study) {
    ggplot(data, aes(x = con, y = .data[[y_est]], fill = con, color = con)) +
      geom_rain(alpha = .5, rain.side = 'l',
                boxplot.args = list(color = "black", outlier.shape = NA),
                boxplot.args.pos = list(
                  position = ggpp::position_dodgenudge(x = .1), width = 0.1
                )) +
      facet_grid(~model) +
      ylab(y_lab) +
      ylim(y_min,y_max)+
      theme_classic() +
      scale_fill_manual(values = pal) +
      scale_color_manual(values = pal) +
      guides(fill = 'none', color = 'none') +
      theme(text = element_text(family = "Times New Roman"),
            axis.text = element_text(size = 12, angle = 45, hjust = 1),
            axis.title = element_text(size = 12),
            legend.text = element_text(size = 12),
            legend.title = element_text(size = 12),
            plot.title = element_text(size = 16)) +
      ggtitle(study) 
  }
  
  # Create plots for each study using map
  plot_list <- purrr::map(c("mls", "ahrb", "abcd"), ~ dataset %>%
                            filter(study == .x) %>%
                            create_plot(.x))
  
  # Combine plots into a single patchwork object
  combined_plots <- wrap_plots(plotlist = plot_list)
  
  # Add patchwork annotations
  combined_plots 
}

update_mls_fwhm <- function(df) {
  df <- df %>%
    mutate(fwhm = case_when(
      fwhm == "fwhm-3.0" ~ "fwhm-3.6",
      fwhm == "fwhm-4.0" ~ "fwhm-4.8",
      fwhm == "fwhm-5.0" ~ "fwhm-6.0",
      fwhm == "fwhm-6.0" ~ "fwhm-7.2",
      fwhm == "fwhm-7.0" ~ "fwhm-8.4",
      TRUE ~ as.character(fwhm)  # if none of the conditions match, keep the original value
    ))
  return(df)
}

create_specr_plot <- function(summary_df, est_label) {
  plot_a = plot_curve(df = summary_df, ci = TRUE, desc = FALSE, legend = FALSE, null = 0)
  plot_b <- plot_choices(df = summary_df, choices = c("fwhm", "motion","contrast","model"), desc = F, null = 0) +
    labs(y = "Variables", x = paste("Ordered Specification Curve \n",est_label, "coefficient"))
  
  cowplot::plot_grid(plot_a, plot_b, ncol = 1, align = "v", axis = 'tblr',
                     labels = c('A', 'B'), rel_heights = c(1, 2),
                     label_fontfamily = "Times", label_size = 12)
}

calculate_summary_stats <- function(data, variable, est_type) {
  data %>%
    select(study, {{variable}}) %>% 
    group_by(study) %>% 
    summarise(
      "est_type" = est_type,
      "median" = median({{variable}}, na.rm = TRUE),
      "mean" = mean({{variable}}, na.rm = TRUE),
      "sd" = sd({{variable}}, na.rm = TRUE),
      "min" = min({{variable}}, na.rm = TRUE),
      "max" = max({{variable}}, na.rm = TRUE) 
    ) %>% 
  kbl(format = "html", 
      booktabs = TRUE) %>% 
  kable_styling(full_width = FALSE)
  
}

The packages are automatically loaded using pacman. The reported .html was last ran on the system: x86_64-apple-darwin17.0 and R version: R version 4.2.1 (2022-06-23) In the Stage 1 PCI Registered Report we are focused on Individual Continouous (intraclass correlation) and the binary/continuous group similarity (jaccard and spearman). This descriptive file includes the output information for the distribution plots, min/max/mean/median of the estimates and the spec plots for each estimate time. This report is for the between-session estimates

continuous

We stated:

Aim1: the range and distribution of median ICCs across each study (three) and analytic decision category (four) are plotted across suprathreshold task-positive and subthreshold ICCs using Rainclouds (Allen et al., 2019) and the median and standard deviation is reported in a table.

to visualize the ordered median ICCs across the 360 model permutations for suprathreshold task-positive and subthreshold ICCs, specification curve analyses are used (Simonsohn et al., 2020). Specifically, results across the 360 model permutations are reported using a specification curve to represent the range of estimated effects across the variable permutations. This consists of two panels: Panel A represents the ordered ICC coefficients and the ICC’s associated 95% confidence interval colored based on no significance (gray), negative (red) or positive (blue) significance from the Null (Null here is 0) and Panel B represents the analytic decisions from each of the four categories (see Table 1) that produced the median ICC estimates. In the main text, to compare the highest and lowest ICC’s produced by the model permutations, the 25th percentile and 75th percentile median ICC estimates from the 360 models are reported separately for suprathreshold task-positive and subthreshold activation (the specification curve for all ICC estimates for suprathreshold task-positive and subthreshold activation are provided as supplemental information).

Aim2: the range and distribution of median Between Subject and Within Subject across each study and analytic decision category are plotted across suprathreshold task-positive and subthreshold ICCs using Rainclouds.

two separate specification curve analyses report the ordered median Between Subject and Within Subject coefficients in one panel and the analytic decisions that produced the Between Subject and Within Subject estimates in a second panel separately for suprathreshold task-positive and subthreshold activation

group similarity

We stated:

Aim1: For each study, the coefficients are plotted to reflect the distribution and range of coefficients. Both Jaccard’s and Spearman correlation are reported separately.

1 Load data

1.1 ABCD

1.2 AHRB

1.3 MLS

1.4 combine data

2 Sample descriptives

2.1 ABCD

# ABCD NDA DATA Info
abcd_nda$subject <- abcd_nda$participant_id
abcd_nda$subject <- gsub("_", "", abcd_nda$subject)
abcd_nda$session <- abcd_nda$eventname
abcd_nda$session <- gsub('2_year_follow_up_y_arm_1', '2YearFollowUpYArm1', abcd_nda$session)
abcd_nda$session <- gsub('baseline_year_1_arm_1', 'baselineYear1Arm1', abcd_nda$session)

# subsat subjects that match final QC'd listed
abcd_nda_subset <- abcd_nda %>% 
  filter(subject %in% abcd_site13_ids$V1) %>% 
  select(!c(subjectkey,eventname))

abcd_beh$subject <- gsub("sub-", "", abcd_beh$subject)
abcd_beh_r <- abcd_beh %>% 
  filter(subject %in% abcd_site13_ids$V1)
abcd_beh_r$sample <- "ABCD"

# estimate days
abcd_days <- abcd_nda_subset %>% 
  mutate(date_r = as.Date(interview_date, format = "%m/%d/%Y")) %>% 
  group_by(subject) %>% 
  arrange(date_r) %>% 
  mutate(days_btwn_scans =  as.integer(diff(date_r))) %>% 
  select(days_btwn_scans, subject, interview_age,sex, race_ethnicity) 

colnames(abcd_days) <- c("days","subject","Age","Sex","Race")
abcd_days$sample <- "ABCD"

abcd_days <- abcd_days %>% 
  mutate(Race = case_when(
      Race == 1 ~ "White",
      Race == 2 ~ "Black",
      Race == 3 ~ "Hispanic",
      Race == 4 ~ "Asian",
      Race == 5 ~ "Other",
      TRUE ~ as.character(Race)  # if none of the conditions match, keep the original value
    ))

rm(abcd_nda,abcd_site13_ids,abcd_beh)

2.2 AHRB

ahrb_beh_r <- ahrb_beh %>% 
  filter(subject %in% ahrb_ids$V1)
ahrb_beh_r$sample <- "AHRB"

ahrb_demo_subset <- ahrb_demo %>% 
  filter(participant_id %in% ahrb_ids$V1) %>% 
  mutate(days = days_ses1toses2)
ahrb_demo_subset$sample <- "AHRB"

ahrb_demo_subset <- ahrb_demo_subset %>% 
  mutate(Race = case_when(
      race == 1 ~ "White",
      race == 2 ~ "Black",
      race == 3 ~ "Hispanic",
      race == 4 ~ "Other",
      TRUE ~ as.character(race)  # if none of the conditions match, keep the original value
    ))

rm(ahrb_demo, ahrb_ids,ahrb_beh)

2.3 mls

mls_beh_r <- mls_beh %>% 
  filter(subject %in% mls_ids$V1)
mls_beh_r$sample <- "MLS"

mls_demo_subset <- mls_demo %>% 
  filter(participant_id %in% mls_ids$V1) %>% 
  mutate(days = if_else(is.na(reward_days_ses1toses2),
                        reward_ses2toses3,reward_days_ses1toses2)
  ) %>% 
  filter(session==1) %>% 
  rename()
mls_demo_subset$sample <- "MLS"
mls_demo_subset <- mls_demo_subset %>% 
  mutate(Race = case_when(
      Race == 1 ~ "White",
      Race == 2 ~ "Black",
      Race == 13 ~ "Hispanic",
      TRUE ~ "Other"  # if none of the conditions match, keep the original value
    ))

rm(mls_demo, mls_ids, mls_beh)
dur_scans <- rbind(abcd_days %>% select(subject,days,sample) %>% unique(),
                   ahrb_demo_subset %>% rename("subject" = participant_id) %>% select(subject,days,sample),
                   mls_demo_subset %>% rename("subject" = participant_id) %>% select(subject,days,sample))

2.4 Days ses-to-ses

Plot distribution of days between scans

dur_scans %>% 
  ggplot(aes(x = subject, y = days, fill = sample, color="white")) +
  geom_bar(stat = 'identity', position = 'dodge') +
  ylab("Days between scans") +
  xlab("") +
  scale_fill_manual(values = pal) +
  scale_color_manual(values = "white") +
  theme_minimal() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), text = element_text(family = "Times New Roman")) +
  facet_grid(~sample, scales = 'free_x') +
  guides(fill = FALSE, color = FALSE, scale = "none")
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

dur_scans %>%
    group_by(sample) %>% 
    summarise(
      "Sample (N)" = n(),
      "Avg Days Between Scans" = mean(days),
      "SD Days Between Scans" = sd(days)
    ) %>% 
  kbl(format = "html", 
      booktabs = TRUE) %>% 
  kable_styling(full_width = FALSE)
sample Sample (N) Avg Days Between Scans SD Days Between Scans
ABCD 119 747.1597 79.09381
AHRB 60 419.4167 80.13660
MLS 81 1088.7037 623.50815

2.5 demographics

demo_comb <- rbind(abcd_days %>% distinct(subject, .keep_all = TRUE) %>% 
                     mutate(sample = "ABCD", Age = round(Age/12,2)) %>% 
                     select(sample, subject,Age, Sex, Race, days),
                   ahrb_demo_subset %>% 
                     mutate(sample = "AHRB", Sex = if_else(sex==0,"F","M")) %>% 
                     rename("Age"= age, "subject" = participant_id) %>% 
                     select(sample, subject,Age, Sex, Race, days),
                   mls_demo_subset %>% 
                     mutate(sample = "MLS", Age = round(if_else(is.na(reward_days_ses1toses2),
                                                                ((ScanAge * 365) + reward_ses2toses3) / 365,
                                                                ScanAge),1),
                            Sex = if_else(Sex==1,"F","M")) %>% 
                     rename("subject" = participant_id) %>% 
                     select(sample, subject,Age,Sex, Race, days) 
                   ) 


table1(~Age + factor(Sex) + factor(Race) + days | sample, data = demo_comb)
ABCD
(N=119)		 AHRB
(N=60)		 MLS
(N=81)		 Overall
(N=260)
Age
Mean (SD) 9.79 (0.587) 19.3 (1.30) 20.7 (2.26) 15.4 (5.35)
Median [Min, Max] 9.83 [9.00, 11.0] 19.2 [17.2, 21.4] 20.2 [18.0, 26.8] 17.8 [9.00, 26.8]
factor(Sex)
F 58 (48.7%) 35 (58.3%) 31 (38.3%) 124 (47.7%)
M 61 (51.3%) 25 (41.7%) 50 (61.7%) 136 (52.3%)
factor(Race)
Asian 4 (3.4%) 0 (0%) 0 (0%) 4 (1.5%)
Black 14 (11.8%) 10 (16.7%) 2 (2.5%) 26 (10.0%)
Hispanic 8 (6.7%) 3 (5.0%) 5 (6.2%) 16 (6.2%)
Other 15 (12.6%) 5 (8.3%) 1 (1.2%) 21 (8.1%)
White 78 (65.5%) 42 (70.0%) 73 (90.1%) 193 (74.2%)
days
Mean (SD) 747 (79.1) 419 (80.1) 1090 (624) 778 (430)
Median [Min, Max] 725 [648, 1070] 391 [352, 692] 756 [332, 2970] 712 [332, 2970]

2.6 behavioral descriptives

beh_dat <- rbind(abcd_beh_r,ahrb_beh_r,mls_beh_r) %>% 
   mutate(session = gsub("baselineYear1Arm1", 1, session), session = gsub("2YearFollowUpYArm1", 2, session))
rm(abcd_beh_r,ahrb_beh_r,mls_beh_r)
beh_dat <- beh_dat %>% 
  mutate(mFD = (mFD_r1+mFD_r2)/2,
         avg_acc = (acc_r1+acc_r2)/2,
         avg_mrt = (mrt_r1+mrt_r2)/2)

# mean/sd, min, max
calculate_summary_stats(beh_dat %>% rename("study" = sample) %>% filter(session=="1"), mFD, "Ses-1: Mean FD")
study est_type median mean sd min max
ABCD Ses-1: Mean FD 0.2002500 0.2465254 0.1531422 0.0640000 0.8615000
AHRB Ses-1: Mean FD 0.1169541 0.1249398 0.0430753 0.0537253 0.2428102
MLS Ses-1: Mean FD 0.1000000 0.1026543 0.0340228 0.0400000 0.2450000 
calculate_summary_stats(beh_dat %>% rename("study" = sample) %>% filter(session=="2"), mFD, "Ses-2: Mean FD")
study est_type median mean sd min max
ABCD Ses-2: Mean FD 0.1825000 0.2390254 0.2061568 0.0450000 1.2930000
AHRB Ses-2: Mean FD 0.1194774 0.1429046 0.0769877 0.0596695 0.5073195
MLS Ses-2: Mean FD 0.0850000 0.0932716 0.0334866 0.0500000 0.2050000 
calculate_summary_stats(beh_dat %>% rename("study" = sample) %>% filter(session=="1"), avg_acc, 
                        "Ses-1: Avg Probe Acc %")
study est_type median mean sd min max
ABCD Ses-1: Avg Probe Acc % 0.55 0.5484746 0.0418101 0.44 0.63
AHRB Ses-1: Avg Probe Acc % 0.58 0.5718333 0.0370749 0.48 0.66
MLS Ses-1: Avg Probe Acc % 0.72 0.7166049 0.1309457 0.40 0.94 
calculate_summary_stats(beh_dat %>% rename("study" = sample) %>% filter(session=="2"), avg_acc, 
                        "Ses-1: Avg Probe Acc %")
study est_type median mean sd min max
ABCD Ses-1: Avg Probe Acc % 0.57 0.5636441 0.0368333 0.440 0.63
AHRB Ses-1: Avg Probe Acc % 0.59 0.5808333 0.0302695 0.510 0.65
MLS Ses-1: Avg Probe Acc % 0.68 0.6685802 0.1208699 0.365 0.94 
calculate_summary_stats(beh_dat %>% rename("study" = sample) %>% filter(session=="1"), avg_mrt, 
                        "Ses-1: Avg Probe MRT (ms)")
study est_type median mean sd min max
ABCD Ses-1: Avg Probe MRT (ms) 303.4650 306.7928 34.49412 233.465 406.2150
AHRB Ses-1: Avg Probe MRT (ms) 296.6986 297.1392 18.47109 236.847 337.7957
MLS Ses-1: Avg Probe MRT (ms) 200.1000 204.3228 28.93196 146.835 268.1850 
calculate_summary_stats(beh_dat %>% rename("study" = sample) %>% filter(session=="2"), avg_mrt, 
                        "Ses-2: Avg Probe MRT (ms)")
study est_type median mean sd min max
ABCD Ses-2: Avg Probe MRT (ms) 270.7000 274.3485 34.13183 214.2600 439.295
AHRB Ses-2: Avg Probe MRT (ms) 243.0866 248.4726 21.47633 217.6078 313.270
MLS Ses-2: Avg Probe MRT (ms) 215.5700 210.0565 30.03950 105.7750 277.430 
acc_plt <- beh_dat %>% 
  ggplot(aes(x = sample, y = avg_acc, fill = sample, color = sample)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .1), width = 0.1
            )) +
  facet_grid(~session) +
  theme_classic() +
  labs(#title = paste("Distribution of Probe Accuracy (%) Across Sample and Sessions"),
          x = "Sample", y = "Probe Acc (%)") +
  scale_fill_manual(values = pal) +
  scale_color_manual(values = pal) +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 12, angle = 45, hjust = 1),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        plot.title = element_text(size = 16))

mrt_plt <- beh_dat %>% 
  ggplot(aes(x = sample, y = avg_mrt, fill = sample, color = sample)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .1), width = 0.1
            )) +
  facet_grid(~session) +
  theme_classic() +
  labs(#title = "Distribution of Probe Mean RT (ms) Across Sample and Sessions",
          x = "Sample", y = "MRT (ms)") +
  scale_fill_manual(values = pal) +
  scale_color_manual(values = pal) +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 12, angle = 45, hjust = 1),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        plot.title = element_text(size = 16))

mfd_plt <- beh_dat %>% 
  ggplot(aes(x = sample, y = mFD, fill = sample, color = sample)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .1), width = 0.1
            )) +
  facet_grid(~session) +
  theme_classic() +
  labs(#title = "Distribution of mean Framewise-displace (mFD) Across Sample and Sessions",
          x = "Sample", y = "mFD") +
  ylim(0,1) +
  scale_fill_manual(values = pal) +
  scale_color_manual(values = pal) +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 12, angle = 45, hjust = 1),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        plot.title = element_text(size = 16))

cat("Plotting A = Avg Acc; B = Avg MRT; C = MeanFD ")
## Plotting A = Avg Acc; B = Avg MRT; C = MeanFD
(acc_plt | mrt_plt | mfd_plt) + plot_annotation(tag_levels = c("A","B","C")) & theme(plot.tag = element_text(size = 30, face = "bold"))

3 Plot distributions

Below, running the steps to summarize the different Intraclass correlation (ICC), Mean Squared Between Subject (Between Subject), Mean Square Within Subject (Within Subject), and jaccard and spearman similarty for the model combinations 360 across samples

3.1 ICC

Plotting overall and for each of [four] categories

3.1.1 Subthreshold Mask

Creating rainclouds via ggrain

subset <- by_sample_conmod(icc_subthresh, y_est = "median_est", y_min =-.1, y_max=.6, y_lab = "Median ICC")
fwhm_rg = create_plot(icc_subthresh, y_est = "median_est", type = "fwhm",  thresh = "sub-threshold", 
                      y_min=-.1, y_max=.6, y_lab = "Median ICC", add_title=FALSE)
motion_rg = create_plot(icc_subthresh, y_est = "median_est", type = "motion", thresh = "sub-threshold", 
                        y_min=-.1, y_max=.6, y_lab = "Median ICC", add_title=FALSE)
modeltype_rg = create_plot(icc_subthresh, y_est = "median_est", type = "model", thresh = "sub-threshold",
                           y_min=-.1, y_max=.6, y_lab = "Median ICC", add_title=FALSE)
contrast_rg = create_plot(icc_subthresh, y_est = "median_est", type = "con", thresh = "sub-threshold", 
                          y_min=-.1, y_max=.6, y_lab = "Median ICC", add_title=FALSE)

cat("Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for sub-threshold mask")
## Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for sub-threshold mask
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold")) 

subset

3.1.2 Suprathreshold Mask

subset <- by_sample_conmod(icc_suprathresh, y_est = "median_est", y_min=-.1, y_max=.6, y_lab = "Median ICC")
fwhm_rg = create_plot(icc_suprathresh, y_est = "median_est", type = "fwhm",  thresh = "supra-threshold", 
                      y_min=-.1, y_max=.6, y_lab = "Median ICC", add_title=FALSE)
motion_rg = create_plot(icc_suprathresh, y_est = "median_est", type = "motion", thresh = "supra-threshold", 
                        y_min=-.1, y_max=.6, y_lab = "Median ICC", add_title=FALSE)
modeltype_rg = create_plot(icc_suprathresh, y_est = "median_est", type = "model", thresh = "supra-threshold",
                           y_min=-.1, y_max=.6, y_lab = "Median ICC", add_title=FALSE)
contrast_rg = create_plot(icc_suprathresh, y_est = "median_est", type = "con", thresh = "supra-threshold", 
                          y_min=-.1, y_max=.6, y_lab = "Median ICC", add_title=FALSE)


cat("Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for supra-threshold mask")
## Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for supra-threshold mask
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold")) 

subset

3.1.3 Nacc avg Mask

3.1.3.2 Left

fwhm_rg = create_plot(icc_nacc, y_est = "avg_left", "fwhm","Avg Left NAcc", y_min=-.3, y_max=.6, y_lab = "ICC", add_title=FALSE)
motion_rg = create_plot(icc_nacc, y_est = "avg_left", "motion","Avg Left NAcc", y_min=-.3, y_max=.6, y_lab = "ICC", add_title=FALSE)
modeltype_rg = create_plot(icc_nacc, y_est = "avg_left", "model","Avg Left NAcc",y_min=-.3, y_max=.6, y_lab = "ICC", add_title=FALSE)
contrast_rg = create_plot(icc_nacc, y_est = "avg_left", "con","sAvg Left NAcc", y_min=-.3, y_max=.6, y_lab = "ICC", add_title=FALSE)

cat("Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for Left NAcc")
## Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for Left NAcc
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold"))

3.2 Between Subject

Plotting overall and for each of [four] categories

3.2.1 Subthreshold Mask

subset <- by_sample_conmod(bs_subthresh, y_est = "median_est", y_min=-.01, y_max=2, y_lab = "Median Between Subject")
fwhm_rg = create_plot(bs_subthresh, y_est = "median_est", type = "fwhm",  thresh = "sub-threshold", 
                      y_min=-.01, y_max=2, y_lab = "Median Between Subject", add_title=FALSE)
motion_rg = create_plot(bs_subthresh, y_est = "median_est", type = "motion", thresh = "sub-threshold", 
                        y_min=-.01, y_max=2, y_lab = "Median Between Subject", add_title=FALSE)
modeltype_rg = create_plot(bs_subthresh, y_est = "median_est", type = "model", thresh = "sub-threshold",
                          y_min=-.01, y_max=2, y_lab = "Median Between Subject", add_title=FALSE)
contrast_rg = create_plot(bs_subthresh, y_est = "median_est", type = "con", thresh = "sub-threshold", 
                         y_min=-.01, y_max=2, y_lab = "Median Between Subject", add_title=FALSE)

cat("Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for sub-threshold mask")
## Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for sub-threshold mask
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold")) 

subset

3.2.2 Suprathreshold Mask

subset <- by_sample_conmod(bs_suprathresh, y_est = "median_est", y_min=-1, y_max=2, y_lab = "Median Between Subject")
fwhm_rg = create_plot(bs_suprathresh, y_est = "median_est", type = "fwhm",  thresh = "supra-threshold", 
                      y_min=-.01, y_max=2, y_lab = "Median Between Subject", add_title=FALSE)
motion_rg = create_plot(bs_suprathresh, y_est = "median_est", type = "motion", thresh = "supra-threshold", 
                        y_min=-.01, y_max=2, y_lab = "Median Between Subject", add_title=FALSE)
modeltype_rg = create_plot(bs_suprathresh, y_est = "median_est", type = "model", thresh = "supra-threshold",
                           y_min=-.01, y_max=2, y_lab = "Median Between Subject", add_title=FALSE)
contrast_rg = create_plot(bs_suprathresh, y_est = "median_est", type = "con", thresh = "supra-threshold", 
                          y_min=-.01, y_max=2, y_lab = "Median Between Subject", add_title=FALSE)

cat("Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for supra-threshold mask")
## Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for supra-threshold mask
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold")) 

subset

3.3 Within Subject

Plotting overall and for each of [four] categories

3.3.1 Subthreshold Mask

subset <- by_sample_conmod(ws_subthresh, y_est = "median_est", y_min=0, y_max=7, y_lab = "Median Within Subject")
fwhm_rg = create_plot(ws_subthresh, y_est = "median_est", type = "fwhm",  thresh = "sub-threshold", 
                      y_min=0, y_max=7, y_lab = "Median Within Subject", add_title=FALSE)
motion_rg = create_plot(ws_subthresh, y_est = "median_est", type = "motion", thresh = "sub-threshold", 
                        y_min=0, y_max=7, y_lab = "Median Within Subject", add_title=FALSE)
modeltype_rg = create_plot(ws_subthresh, y_est = "median_est", type = "model", thresh = "sub-threshold",
                           y_min=0, y_max=7, y_lab = "Median Within Subject", add_title=FALSE)
contrast_rg = create_plot(ws_subthresh, y_est = "median_est", type = "con", thresh = "sub-threshold", 
                          y_min=0, y_max=7, y_lab = "Median Within Subject", add_title=FALSE)

cat("Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for sub-threshold mask")
## Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for sub-threshold mask
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold")) 

subset

3.3.2 Suprathreshold Mask

subset <- by_sample_conmod(ws_suprathresh, y_est = "median_est", y_min=0, y_max=7, y_lab = "Median Within Subject")
fwhm_rg = create_plot(ws_suprathresh, y_est = "median_est", type = "fwhm",  thresh = "supra-threshold", 
                      y_min=0, y_max=7, y_lab = "Median Within Subject", add_title=FALSE)
motion_rg = create_plot(ws_suprathresh, y_est = "median_est", type = "motion", thresh = "supra-threshold", 
                        y_min=0, y_max=7, y_lab = "Median Within Subject", add_title=FALSE)
modeltype_rg = create_plot(ws_suprathresh, y_est = "median_est", type = "model", thresh = "supra-threshold",
                           y_min=0, y_max=7, y_lab = "Median Within Subject", add_title=FALSE)
contrast_rg = create_plot(ws_suprathresh, y_est = "median_est", type = "con", thresh = "supra-threshold", 
                          y_min=0, y_max=7, y_lab = "Median Within Subject", add_title=FALSE)

cat("Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for supra-threshold mask")
## Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast for supra-threshold mask
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold")) 

subset

3.4 Similarity

3.4.1 jaccard

subset <- by_sample_conmod(similarity_df, y_est = "jaccard", y_min=0, y_max = 1, y_lab = "Jaccard")

fwhm_rg = create_plot(similarity_df, y_est = "jaccard", type = "fwhm", thresh = "Jaccard", 
                      y_min=0, y_max = 1, y_lab = "Jaccard", add_title=FALSE)
motion_rg = create_plot(similarity_df, y_est = "jaccard", type = "motion", thresh = "Jaccard", 
                        y_min=0, y_max = 1, y_lab = "Jaccard", add_title=FALSE)
modeltype_rg = create_plot(similarity_df, y_est = "jaccard", type = "model", thresh="Jaccrd", 
                           y_min=0, y_max = 1, y_lab = "Jaccard", add_title=FALSE)
contrast_rg = create_plot(similarity_df, y_est = "jaccard", type = "con", thresh = "Jaccard", 
                          y_min=0, y_max = 1, y_lab = "Jaccard", add_title=FALSE)


cat("Jaccard: Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast")
## Jaccard: Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold")) 

subset

3.4.2 spearman

subset <- by_sample_conmod(similarity_df, y_est = "spearman", y_min=0, y_max = 1, y_lab = "Spearman")

fwhm_rg = create_plot(similarity_df, y_est = "spearman", type = "fwhm", thresh = "Spearman", 
                      y_min=0, y_max = 1, y_lab = "Spearman", add_title=FALSE)
motion_rg = create_plot(similarity_df, y_est = "spearman", type = "motion", thresh = "Spearman", 
                        y_min=0, y_max = 1, y_lab = "Spearman", add_title=FALSE)
modeltype_rg = create_plot(similarity_df, y_est = "spearman", type = "model", thresh="Spearman", 
                           y_min=0, y_max = 1, y_lab = "Spearman", add_title=FALSE)
contrast_rg = create_plot(similarity_df, y_est = "spearman", type = "con", thresh = "Spearman", 
                          y_min=0, y_max = 1, y_lab = "Spearman", add_title=FALSE)

cat("Spearman: Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast")
## Spearman: Plotting A = Motion; B = FWHM; C = Model Parameterization, D = Contrast
(motion_rg | fwhm_rg) / (modeltype_rg | contrast_rg) + plot_annotation(tag_levels = c("A","B","C","D")) & theme(plot.tag = element_text(size = 30, face = "bold")) 

subset

3.5 Group Cohen’s ~ ICC

cohens_sim <- gather(similarity_df, key = "Ses", value = "spearman", ses1_icc_cohensd:ses2_icc_cohensd) %>% 
  mutate(Run = case_when(
    Ses == "ses1_icc_cohensd" ~ "ses1",
    Ses == "ses_icc_cohensd" ~ "ses2",
    TRUE ~ Ses
  ))
cat("\tICC, Between Subject and Within Subject median, miniumum and maximum")
##  ICC, Between Subject and Within Subject median, miniumum and maximum
calculate_summary_stats(similarity_df, ses1_icc_cohensd, "ICC ~ Ses1 Cohen's D")
study est_type median mean sd min max
abcd ICC ~ Ses1 Cohen’s D 0.1019650 0.0088166 0.1973574 -0.4001350 0.2888389
ahrb ICC ~ Ses1 Cohen’s D 0.1723025 0.1124787 0.2548283 -0.4450425 0.5265826
mls ICC ~ Ses1 Cohen’s D 0.1213723 0.1168365 0.1887960 -0.2811225 0.4317429 
calculate_summary_stats(similarity_df, ses2_icc_cohensd, "ICC ~ Ses2 Cohen's D")
study est_type median mean sd min max
abcd ICC ~ Ses2 Cohen’s D 0.1060762 -0.0078938 0.2279903 -0.4635707 0.2964393
ahrb ICC ~ Ses2 Cohen’s D 0.2013921 0.1173448 0.2557813 -0.4294658 0.5289596
mls ICC ~ Ses2 Cohen’s D 0.1360693 0.1134261 0.1768921 -0.3098683 0.3929272 
cohens_sim %>% ggplot(aes(x = fwhm, y = spearman, fill = fwhm, color = fwhm)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .2), width = 0.1
            )) +
  facet_grid(~study + Ses) +
  theme_classic() +
  ggtitle("Distribution by FWHM category") +
  scale_fill_manual(values = pal) +
  scale_color_manual(values = pal) +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 12, angle = 45, hjust = 1),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        plot.title = element_text(size = 16))

cohens_sim %>% ggplot(aes(x = con, y = spearman, fill = con, color = con)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .2), width = 0.1
            )) +
  facet_grid(~study + Ses) +
  theme_classic() +
  ggtitle("Distribution by Contrast category") +
  scale_fill_manual(values = pal) +
  scale_color_manual(values = pal) +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 12, angle = 45, hjust = 1),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        plot.title = element_text(size = 16))

cohens_sim %>% ggplot(aes(x = motion, y = spearman, fill = motion, color = motion)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .2), width = 0.1
            )) +
  facet_grid(~study + Ses) +
  theme_classic() +
  ggtitle("Distribution by Motion category") +
  scale_fill_manual(values = pal) +
  scale_color_manual(values = pal) +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 12, angle = 45, hjust = 1),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        plot.title = element_text(size = 16))

cohens_sim %>% ggplot(aes(x = model, y = spearman, fill = model, color = model)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .2), width = 0.1
            )) +
  facet_grid(~study + Ses) +
  theme_classic() +
  ggtitle("Distribution by Model category") +
  scale_fill_manual(values = pal) +
  scale_color_manual(values = pal) +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 12, angle = 45, hjust = 1),
        axis.title = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        plot.title = element_text(size = 16))

3.6 Dev Comp.

3.6.1 Full

older <- icc_suprathresh %>% filter(study %in% c('mls', 'ahrb')) 
younger <- icc_suprathresh %>% filter(study == 'abcd') 

cat("Comparing (x) older samples AHRB/MLS versus (y) younger sample ABCD")
## Comparing (x) older samples AHRB/MLS versus (y) younger sample ABCD
t.test(older$median_est, younger$median_est)
## 
##  Welch Two Sample t-test
## 
## data:  older$median_est and younger$median_est
## t = 11.324, df = 697.29, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.06335638 0.08993529
## sample estimates:
## mean of x mean of y 
## 0.2236833 0.1470375
effectsize::cohens_d(older$median_est, younger$median_est, pooled = TRUE)
## Cohen's d |       95% CI
## ------------------------
## 0.76      | [0.60, 0.92]
## 
## - Estimated using pooled SD.

3.6.2 FWHM

filtered_df <- icc_suprathresh %>%
  filter(motion == 'opt1' & model == 'CueMod' & con == 'LgainNeut')

# Separate the fwhm values based on study groups
older <- filtered_df %>% filter(study %in% c('mls', 'ahrb')) 
younger <- filtered_df %>% filter(study == 'abcd') 

cat("Comparing (x) older samples AHRB/MLS (n = ", nrow(older),")", "versus (y) younger sample ABCD (n = ", nrow(younger),")")
## Comparing (x) older samples AHRB/MLS (n =  10 ) versus (y) younger sample ABCD (n =  5 )
t.test(older$median_est, younger$median_est)
## 
##  Welch Two Sample t-test
## 
## data:  older$median_est and younger$median_est
## t = 0.54854, df = 12.236, p-value = 0.5932
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.04415646  0.07395646
## sample estimates:
## mean of x mean of y 
##    0.1933    0.1784
effectsize::cohens_d(older$median_est, younger$median_est, pooled = TRUE)
## Cohen's d |        95% CI
## -------------------------
## 0.23      | [-0.86, 1.30]
## 
## - Estimated using pooled SD.

3.6.3 Model Parameterization

filtered_df <- icc_suprathresh %>%
  filter(motion == 'opt1' & fwhm == '4.8' & con == 'LgainNeut')

# Separate the fwhm values based on study groups
older <- filtered_df %>% filter(study %in% c('mls', 'ahrb')) 
younger <- filtered_df %>% filter(study == 'abcd') 

cat("Comparing (x) older samples AHRB/MLS (n = ", nrow(older),")", "versus (y) younger sample ABCD (n = ", nrow(younger),")")
## Comparing (x) older samples AHRB/MLS (n =  6 ) versus (y) younger sample ABCD (n =  3 )
t.test(older$median_est, younger$median_est)
## 
##  Welch Two Sample t-test
## 
## data:  older$median_est and younger$median_est
## t = 1.2152, df = 6.8276, p-value = 0.2646
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0331380  0.1024713
## sample estimates:
## mean of x mean of y 
## 0.1850000 0.1503333
effectsize::cohens_d(older$median_est, younger$median_est, pooled = TRUE)
## Cohen's d |        95% CI
## -------------------------
## 0.64      | [-0.80, 2.05]
## 
## - Estimated using pooled SD.

3.6.4 Contrasts

filtered_df <- icc_suprathresh %>%
  filter(motion == 'opt1' & fwhm == '3.6' & model == 'FixMod')

# Separate the fwhm values based on study groups
older <- filtered_df %>% filter(study %in% c('mls', 'ahrb')) 
younger <- filtered_df %>% filter(study == 'abcd') 

cat("Comparing (x) older samples AHRB/MLS (n = ", nrow(older),")", "versus (y) younger sample ABCD (n = ", nrow(younger),")")
## Comparing (x) older samples AHRB/MLS (n =  8 ) versus (y) younger sample ABCD (n =  4 )
t.test(older$median_est, younger$median_est)
## 
##  Welch Two Sample t-test
## 
## data:  older$median_est and younger$median_est
## t = 1.9067, df = 9.9185, p-value = 0.08591
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01522676  0.19447676
## sample estimates:
## mean of x mean of y 
##  0.203375  0.113750
effectsize::cohens_d(older$median_est, younger$median_est, pooled = TRUE)
## Cohen's d |        95% CI
## -------------------------
## 0.90      | [-0.38, 2.14]
## 
## - Estimated using pooled SD.

3.6.5 Motion

filtered_df <- icc_suprathresh %>%
  filter(con == 'LgainNeut' & fwhm == '3.6' & model == 'FixMod')


# Separate the fwhm values based on study groups
older <- filtered_df %>% filter(study %in% c('mls', 'ahrb')) 
younger <- filtered_df %>% filter(study == 'abcd') 

cat("Comparing (x) older samples AHRB/MLS (n = ", nrow(older),")", "versus (y) younger sample ABCD (n = ", nrow(younger),")")
## Comparing (x) older samples AHRB/MLS (n =  8 ) versus (y) younger sample ABCD (n =  4 )
t.test(older$median_est, younger$median_est)
## 
##  Welch Two Sample t-test
## 
## data:  older$median_est and younger$median_est
## t = 2.0015, df = 9.4033, p-value = 0.07501
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.004669914  0.080669914
## sample estimates:
## mean of x mean of y 
##    0.1325    0.0945
effectsize::cohens_d(older$median_est, younger$median_est, pooled = TRUE)
## Cohen's d |        95% CI
## -------------------------
## 1.03      | [-0.27, 2.29]
## 
## - Estimated using pooled SD.

4 Med/Min/Max Across Models

4.1 Subtreshold

cat("Subthreshold mask")
## Subthreshold mask
cat("\tICC, Between Subject and Within Subject median, miniumum and maximum")
##  ICC, Between Subject and Within Subject median, miniumum and maximum
calculate_summary_stats(icc_subthresh, median_est, "ICC")
study est_type median mean sd min max
abcd ICC 0.086 0.0968833 0.0456961 0.024 0.253
ahrb ICC 0.126 0.1526667 0.0878269 0.041 0.402
mls ICC 0.131 0.1435792 0.0665285 0.037 0.320 
calculate_summary_stats(bs_subthresh, median_est, "Between Subject")
study est_type median mean sd min max
abcd Between Subject 0.0335 0.0577125 0.0646393 0.002 0.322
ahrb Between Subject 0.0250 0.0510167 0.0610518 0.003 0.272
mls Between Subject 0.0555 0.1300792 0.1510471 0.002 0.574 
calculate_summary_stats(ws_subthresh, median_est, "Within Subject")
study est_type median mean sd min max
abcd Within Subject 0.4015 0.4821625 0.3360764 0.047 1.877
ahrb Within Subject 0.1895 0.2334250 0.1610097 0.027 0.932
mls Within Subject 0.4655 0.5955125 0.4781952 0.041 2.076

4.2 Suprathreshold

cat("Suprathreshold mask")
## Suprathreshold mask
cat("\tICC, Between Subject and Within Subject median, miniumum and maximum")
##  ICC, Between Subject and Within Subject median, miniumum and maximum
calculate_summary_stats(icc_suprathresh, median_est, "ICC")
study est_type median mean sd min max
abcd ICC 0.1375 0.1470375 0.0674774 0.029 0.318
ahrb ICC 0.2090 0.2310792 0.1277304 0.043 0.532
mls ICC 0.2055 0.2162875 0.0969472 0.058 0.473 
calculate_summary_stats(bs_suprathresh, median_est, "Between Subject")
study est_type median mean sd min max
abcd Between Subject 0.0370 0.0561250 0.0570199 0.002 0.294
ahrb Between Subject 0.0280 0.0526083 0.0606377 0.002 0.278
mls Between Subject 0.0615 0.1352125 0.1572652 0.003 0.683 
calculate_summary_stats(ws_suprathresh, median_est, "Within Subject")
study est_type median mean sd min max
abcd Within Subject 0.2590 0.3002708 0.1996876 0.031 1.097
ahrb Within Subject 0.1225 0.1462917 0.0985820 0.017 0.555
mls Within Subject 0.2790 0.3634750 0.2890710 0.024 1.238

4.3 simlarity

cat("Similarity median, minimum and maximum for jaccard and spearman")
## Similarity median, minimum and maximum for jaccard and spearman
calculate_summary_stats(similarity_df, jaccard, "Jaccard")
study est_type median mean sd min max
abcd Jaccard 0.2488132 0.2571102 0.1301362 0.0208562 0.6046645
ahrb Jaccard 0.3037331 0.3224565 0.1853570 0.0395941 0.7281375
mls Jaccard 0.4178440 0.4347195 0.1234524 0.1979287 0.7410373 
calculate_summary_stats(similarity_df, spearman, "Spearman")
study est_type median mean sd min max
abcd Spearman 0.8029798 0.7568178 0.1322711 0.4012253 0.9350044
ahrb Spearman 0.8243898 0.7428236 0.2085670 0.3161175 0.9695884
mls Spearman 0.8712244 0.8478988 0.0910636 0.5911072 0.9723035

5 Specification Curve

Creating data in a format that is compatible with specr. Needs: estimate (i.e., ICC), std.error, conf.high, conf.low.

5.1 ICC

5.1.1 suprathreshold

5.1.1.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

# first, combine independent model vars into string to create average for each model type
icc_suprathresh$model_type <- paste(icc_suprathresh$fwhm, icc_suprathresh$motion,
                                    icc_suprathresh$con,  icc_suprathresh$model,
                                    sep = "_")

# calculate the avg estimate of ICC across study, standard error and +/- 95% confidence interval. In complete version
df_summ <- icc_suprathresh %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(median_est), std.error = sd(median_est)/sqrt(length(median_est))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","model"), sep = "_|-", remove = FALSE)

# plot
est_label = "supra-threshod ICC"
create_specr_plot(df_summ, est_label)

5.1.1.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subset plots
est_label = "supra-threshold ICC"
create_specr_plot(df_summ_subset, est_label)

5.1.2 subthreshold

5.1.2.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

icc_subthresh$model_type <- paste(icc_subthresh$fwhm, icc_subthresh$motion,
                                  icc_subthresh$con,  icc_subthresh$model,
                                  sep = "_")

df_summ <- icc_subthresh %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(median_est), std.error = sd(median_est)/sqrt(length(median_est))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","model"), sep = "_|-", remove = FALSE)

# plot
est_label = "sub-threshold ICC"
create_specr_plot(df_summ, est_label)

5.1.2.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subset plots
est_label = "sub-threshold ICC"
create_specr_plot(df_summ_subset, est_label)

5.2 Between Subject

5.2.1 suprathreshold

5.2.1.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

bs_suprathresh$model_type <- paste(bs_suprathresh$fwhm, bs_suprathresh$motion,
                                   bs_suprathresh$con,  bs_suprathresh$model,
                                   sep = "_")

df_summ <- bs_suprathresh %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(median_est), std.error = sd(median_est)/sqrt(length(median_est))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","model"), sep = "_|-", remove = FALSE)

# plot
est_label = "supra-threshold Between Subject"
create_specr_plot(df_summ, est_label)

5.2.1.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subset plots
est_label = "supra-threshold Between Subject"
create_specr_plot(df_summ_subset, est_label)

5.2.2 subthreshold

5.2.2.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

bs_subthresh$model_type <- paste(bs_subthresh$fwhm, bs_subthresh$motion,
                                 bs_subthresh$con,  bs_subthresh$model,
                                 sep = "_")


df_summ <- bs_subthresh %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(median_est), std.error = sd(median_est)/sqrt(length(median_est))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","model"), sep = "_|-", remove = FALSE)

# plot
est_label = "sub-threshold Between Subject"
create_specr_plot(df_summ, est_label)

5.2.2.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subset plots
est_label = "suprathreshold Between Subject"
create_specr_plot(df_summ_subset, est_label)

5.3 Within Subject

5.3.1 suprathreshold

5.3.1.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

ws_suprathresh$model_type <- paste(ws_suprathresh$fwhm, ws_suprathresh$motion,
                                   bs_suprathresh$con,  ws_suprathresh$model,
                                   sep = "_")

df_summ <- ws_suprathresh %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(median_est), std.error = sd(median_est)/sqrt(length(median_est))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","model"), sep = "_|-", remove = FALSE)

# plot
est_label = "supra-threshold Within Subject"
create_specr_plot(df_summ, est_label)

5.3.1.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subset plots
est_label = "supra-threshold Within Subject"
create_specr_plot(df_summ_subset, est_label)

5.3.2 subthreshold

5.3.2.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

ws_subthresh$model_type <- paste(ws_subthresh$fwhm, ws_subthresh$motion,
                                 ws_subthresh$con,  ws_subthresh$model,
                                 sep = "_")

# calculate the avg estimate of ICC across study, standard error and +/- 95% confidence interval. In complete version
df_summ <- ws_subthresh %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(median_est), std.error = sd(median_est)/sqrt(length(median_est))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","model"), sep = "_|-", remove = FALSE)

# plot
est_label = "sub-threshold Within Subject"
create_specr_plot(df_summ, est_label)

5.3.2.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subplot
est_label = "sub-threshold Within Subject"
create_specr_plot(df_summ_subset, est_label)

5.4 Similarity

5.4.1 Jaccard

5.4.1.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

similarity_df$model_type <- paste(similarity_df$fwhm, similarity_df$motion,
                                  similarity_df$con, similarity_df$model,
                                  sep = "_")

df_summ <- similarity_df %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(jaccard), std.error = sd(jaccard)/sqrt(length(jaccard))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","model"), sep = "_|-", remove = FALSE)

# plot
est_label = "Jaccard"
create_specr_plot(df_summ, est_label)

5.4.1.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subset plots
est_label = "Jaccard"
create_specr_plot(df_summ_subset, est_label)

5.4.2 Spearman

5.4.2.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

similarity_df$model_type <- paste(similarity_df$fwhm, similarity_df$motion,
                                  similarity_df$con,  similarity_df$model, 
                                  sep = "_")

# calculate the avg estimate of ICC across study, standard error and +/- 95% confidence interval. In complete version
df_summ <- similarity_df %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(spearman), std.error = sd(spearman)/sqrt(length(spearman))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","model"), sep = "_|-", remove = FALSE)

# plot
est_label = "Spearman"
create_specr_plot(df_summ, est_label)

5.4.2.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subset plots
est_label = "Spearman"
create_specr_plot(df_summ_subset, est_label)